/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.io.Serializable;

import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.ThinPlateElement;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;

/** 
 * A class mounts the stiffness matrix of one isotropic element CKZ
 * (triagular plate-bending element nonconforming of 
 * three node developed by Cheung, king and Zienkiewicz), 
 * according to: Cheung, Y. K., king, I. P. and Zienkiewicz, O. C. -
 * "Slab bridges with arbitrary shape and support conditions:
 * a general method of analysis based on finit elements", Proc. Inst.
 * Civil Engng., Vol. 40, pp. 9-36, 1968.
 *
 * @author Samir
 * @since 15/09/2006
 */
public class CKZProblemDriver extends KirchhoffThinPlate implements Serializable {

	private static final long serialVersionUID = 1L;

	/* (non-Javadoc)
	 * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
	 */
	public IMatrix getC(Element e) {
	    int numberOfDegreesOfFreedom = ((ThinPlateElement)e).getNumberOfDegreesOfFreedom();
		IMatrix finalStiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    finalStiffnessMatrix.setZero();
	    IMatrix cartesianNodalCoordsMatrix = ((ThinPlateElement)e).getCartesianNodalCoordsMatrix();
	    IVector centroid = this.getCentroid(cartesianNodalCoordsMatrix);
	    
	    double x1 = cartesianNodalCoordsMatrix.getElement(0,0) - centroid.getElement(0);
	    double x2 = cartesianNodalCoordsMatrix.getElement(1,0) - centroid.getElement(0);
	    double x3 = cartesianNodalCoordsMatrix.getElement(2,0) - centroid.getElement(0);

	    double y1 = cartesianNodalCoordsMatrix.getElement(0,1) - centroid.getElement(1);
	    double y2 = cartesianNodalCoordsMatrix.getElement(1,1) - centroid.getElement(1);
	    double y3 = cartesianNodalCoordsMatrix.getElement(2,1) - centroid.getElement(1);
	    	    
	    double t = ((ThinPlateElement)e).getThickness();
	    double E = ((ThinPlateElement)e).getYoungsModulus();
	    double v = ((ThinPlateElement)e).getPoissonsRatio();
	    
	    IMatrix B = this.getB(x1,x2,x3,y1,y2,y3);
	    IMatrix T = this.getT(x1,x2,x3,y1,y2,y3);
	    IMatrix ADA = this.getADA(x1,x2,x3,y1,y2,y3,E,v,t);
	    
		double A = 0.5*(x2*y3+x3*y1+x1*y2-x2*y1-x1*y3-x3*y2);
	    double cte = 1/(64*Math.pow(A,6));
	    
	    IMatrix BT = new IMatrix(B.getNumRow(),T.getNumCol());
	    BT.setZero();
	    BT.mul(B,T);
	    
	    IMatrix K1 = new IMatrix(BT.getNumCol(),ADA.getNumCol());
	    K1.mulTransposeLeft(BT,ADA);
	    
	    finalStiffnessMatrix.mul(K1,BT);
	    finalStiffnessMatrix.setScale(cte);
	    
	    return(this.signal(this.adjustMatrix(finalStiffnessMatrix)));
	}
	
	/**
	 * Mounts the geometric properties matrix.
	 * @param x1 nodal coordinate of the corner 1 in direction x in cartesian system. 
	 * @param x2 nodal coordinate of the corner 2 in direction x in cartesian system. 
	 * @param x3 nodal coordinate of the corner 3 in direction x in cartesian system. 
	 * @param y1 nodal coordinate of the corner 1 in direction y in cartesian system. 
	 * @param y2 nodal coordinate of the corner 2 in direction y in cartesian system. 
	 * @param y3 nodal coordinate of the corner 3 in direction y in cartesian system. 
	 * @return B matrix.
	 */
	private IMatrix getB(double x1, double x2, double x3, double y1, double y2, double y3){
		double b1 = y2 - y3;   double b2 = y3 - y1;   double b3 = y1 - y2;
		double c1 = x3 - x2;   double c2 = x1 - x3;   double c3 = x2 - x1;
		
		double A = 0.5*(x2*y3+x3*y1+x1*y2-x2*y1-x1*y3-x3*y2);
		double a = A/1.5;
		
		double X11 = b1*b1*a + 2*a*b1*b2;
		double X21 = 2*b1*c1*a + 2*b2*c1*a + 2*b1*c2*a;
		double X31 = c1*c1*a + 2*a*c1*c2;
		double X41 = b1*b1*b2;
		double X51 = 2*b1*c1*b2 + b1*b1*c2;
		double X61 = c1*c1*b2 + 2*b1*c1*c2;
		double X71 = c1*c1*c2;
		
		double X12 = b2*b2*a + 2*a*b2*b3;
		double X22 = 2*b2*c2*a + 2*b3*c2*a + 2*b2*c3*a;
		double X32 = c2*c2*a + 2*a*c2*c3;
		double X42 = b2*b2*b3;
		double X52 = 2*b2*c2*b3 + b2*b2*c3;
		double X62 = c2*c2*b3 + 2*b2*c2*c3;
		double X72 = c2*c2*c3;
		
		double X13 = b3*b3*a + 2*a*b3*b1;
		double X23 = 2*b3*c3*a + 2*b1*c3*a + 2*b3*c1*a;
		double X33 = c3*c3*a + 2*a*c3*c1;
		double X43 = b3*b3*b1;
		double X53 = 2*b3*c3*b1 + b3*b3*c1;
		double X63 = c3*c3*b1 + 2*b3*c3*c1;
		double X73 = c3*c3*c1;
		
		double Y11 = b1*b1*a + 2*a*b1*b3;
		double Y21 = 2*b1*c1*a + 2*b3*c1*a + 2*b1*c3*a;
		double Y31 = c1*c1*a + 2*a*c1*c3;
		double Y41 = b1*b1*b3;
		double Y51 = 2*b1*c1*b3 + b1*b1*c3;
		double Y61 = c1*c1*b3 + 2*b1*c1*c3;
		double Y71 = c1*c1*c3;

		double Y12 = b2*b2*a + 2*a*b2*b1;
		double Y22 = 2*b2*c2*a + 2*b1*c2*a + 2*b2*c1*a;
		double Y32 = c2*c2*a + 2*a*c2*c1;
		double Y42 = b2*b2*b1;
		double Y52 = 2*b2*c2*b1 + b2*b2*c1;
		double Y62 = c2*c2*b1 + 2*b2*c2*c1;
		double Y72 = c2*c2*c1;

		double Y13 = b3*b3*a + 2*a*b3*b2;
		double Y23 = 2*b3*c3*a + 2*b2*c3*a + 2*b3*c2*a;
		double Y33 = c3*c3*a + 2*a*c3*c2;
		double Y43 = b3*b3*b2;
		double Y53 = 2*b3*c3*b2 + b3*b3*c2;
		double Y63 = c3*c3*b2 + 2*b3*c3*c2;
		double Y73 = c3*c3*c2;
		
		double E1 = a*(b1*b2 + b2*b3 + b3*b1);
		double E2 = a*((c1*b2 + b1*c2) + (c2*b3 + b2*c3) + (c3*b1 + b3*c1));
		double E3 = a*(c1*c2 + c2*c3 + c3*c1);
		double E4 = b1*b2*b3;
		double E5 = c1*b2*b3 + c2*b3*b1 + c3*b1*b2;
		double E6 = c1*c2*b3 + c2*c3*b1 + c3*c1*b2;
		double E7 = c1*c2*c3;
		
		double F1 = 0.5*(b3-b2);
		double G1 = 0.5*(c3-c2);
		
		double F2 = 0.5*(b1-b3);
		double G2 = 0.5*(c1-c3);
		
		double F3 = 0.5*(b2-b1);
		double G3 = 0.5*(c2-c1);
		
		double BX11 = X11*b3 - Y11*b2 + E1*F1;
		double BX21 = X21*b3 - Y21*b2 + E2*F1;
		double BX31 = X31*b3 - Y31*b2 + E3*F1;
		double BX41 = X41*b3 - Y41*b2 + E4*F1;
		double BX51 = X51*b3 - Y51*b2 + E5*F1;
		double BX61 = X61*b3 - Y61*b2 + E6*F1;
		double BX71 = X71*b3 - Y71*b2 + E7*F1;
		
		double BX12 = X12*b1 - Y12*b3 + E1*F2;
		double BX22 = X22*b1 - Y22*b3 + E2*F2;
		double BX32 = X32*b1 - Y32*b3 + E3*F2;
		double BX42 = X42*b1 - Y42*b3 + E4*F2;
		double BX52 = X52*b1 - Y52*b3 + E5*F2;
		double BX62 = X62*b1 - Y62*b3 + E6*F2;
		double BX72 = X72*b1 - Y72*b3 + E7*F2;
		
		double BX13 = X13*b2 - Y13*b1 + E1*F3;
		double BX23 = X23*b2 - Y23*b1 + E2*F3;
		double BX33 = X33*b2 - Y33*b1 + E3*F3;
		double BX43 = X43*b2 - Y43*b1 + E4*F3;
		double BX53 = X53*b2 - Y53*b1 + E5*F3;
		double BX63 = X63*b2 - Y63*b1 + E6*F3;
		double BX73 = X73*b2 - Y73*b1 + E7*F3;
		
		double BY11 = X11*c3 - Y11*c2 + E1*G1;
		double BY21 = X21*c3 - Y21*c2 + E2*G1;
		double BY31 = X31*c3 - Y31*c2 + E3*G1;
		double BY41 = X41*c3 - Y41*c2 + E4*G1;
		double BY51 = X51*c3 - Y51*c2 + E5*G1;
		double BY61 = X61*c3 - Y61*c2 + E6*G1;
		double BY71 = X71*c3 - Y71*c2 + E7*G1;
		
		double BY12 = X12*c1 - Y12*c3 + E1*G2;
		double BY22 = X22*c1 - Y22*c3 + E2*G2;
		double BY32 = X32*c1 - Y32*c3 + E3*G2;
		double BY42 = X42*c1 - Y42*c3 + E4*G2;
		double BY52 = X52*c1 - Y52*c3 + E5*G2;
		double BY62 = X62*c1 - Y62*c3 + E6*G2;
		double BY72 = X72*c1 - Y72*c3 + E7*G2;
		
		double BY13 = X13*c2 - Y13*c1 + E1*G3;
		double BY23 = X23*c2 - Y23*c1 + E2*G3;
		double BY33 = X33*c2 - Y33*c1 + E3*G3;
		double BY43 = X43*c2 - Y43*c1 + E4*G3;
		double BY53 = X53*c2 - Y53*c1 + E5*G3;
		double BY63 = X63*c2 - Y63*c1 + E6*G3;
		double BY73 = X73*c2 - Y73*c1 + E7*G3;
		
		IMatrix B = new IMatrix(7,6);
		B.setZero();
		
		B.setElement(0,0,BX11); B.setElement(0,1,BY11); B.setElement(0,2,BX12); B.setElement(0,3,BY12); B.setElement(0,4,BX13); B.setElement(0,5,BY13);
		B.setElement(1,0,BX21); B.setElement(1,1,BY21); B.setElement(1,2,BX22); B.setElement(1,3,BY22); B.setElement(1,4,BX23); B.setElement(1,5,BY23);
		B.setElement(2,0,BX31); B.setElement(2,1,BY31); B.setElement(2,2,BX32); B.setElement(2,3,BY32); B.setElement(2,4,BX33); B.setElement(2,5,BY33);
		B.setElement(3,0,BX41); B.setElement(3,1,BY41); B.setElement(3,2,BX42); B.setElement(3,3,BY42); B.setElement(3,4,BX43); B.setElement(3,5,BY43);
		B.setElement(4,0,BX51); B.setElement(4,1,BY51); B.setElement(4,2,BX52); B.setElement(4,3,BY52); B.setElement(4,4,BX53); B.setElement(4,5,BY53);
		B.setElement(5,0,BX61); B.setElement(5,1,BY61); B.setElement(5,2,BX62); B.setElement(5,3,BY62); B.setElement(5,4,BX63); B.setElement(5,5,BY63);
		B.setElement(6,0,BX71); B.setElement(6,1,BY71); B.setElement(6,2,BX72); B.setElement(6,3,BY72); B.setElement(6,4,BX73); B.setElement(6,5,BY73);
		
		return B;
	}
	
	/**
	 * Mounts the transformation matrix.
	 * @param x1 nodal coordinate of the corner 1 in direction x in cartesian system. 
	 * @param x2 nodal coordinate of the corner 2 in direction x in cartesian system. 
	 * @param x3 nodal coordinate of the corner 3 in direction x in cartesian system. 
	 * @param y1 nodal coordinate of the corner 1 in direction y in cartesian system. 
	 * @param y2 nodal coordinate of the corner 2 in direction y in cartesian system. 
	 * @param y3 nodal coordinate of the corner 3 in direction y in cartesian system. 
	 * @return T matrix.
	 */
	private IMatrix getT(double x1, double x2, double x3, double y1, double y2, double y3){
		double b1 = y2 - y3;   double b2 = y3 - y1;   double b3 = y1 - y2;
		double c1 = x3 - x2;   double c2 = x1 - x3;   double c3 = x2 - x1;
		
		double A = 0.5*(x2*y3+x3*y1+x1*y2-x2*y1-x1*y3-x3*y2);
		
		IMatrix T = new IMatrix(6,9);
		T.setZero();
		
		T.setElement(0,0,(c1/(2*A)));  T.setElement(0,1,1); T.setElement(0,3,(c2/(2*A)));  T.setElement(0,6,(c3/(2*A)));
		T.setElement(1,0,(-b1/(2*A))); T.setElement(1,2,1); T.setElement(1,3,(-b2/(2*A))); T.setElement(1,6,(-b3/(2*A)));
		
		T.setElement(2,0,(c1/(2*A)));  T.setElement(2,3,(c2/(2*A)));  T.setElement(2,4,1); T.setElement(2,6,(c3/(2*A)));
		T.setElement(3,0,(-b1/(2*A))); T.setElement(3,3,(-b2/(2*A))); T.setElement(3,5,1); T.setElement(3,6,(-b3/(2*A)));
		
		T.setElement(4,0,(c1/(2*A)));  T.setElement(4,3,(c2/(2*A)));  T.setElement(4,6,(c3/(2*A)));  T.setElement(4,7,1); 
		T.setElement(5,0,(-b1/(2*A))); T.setElement(5,3,(-b2/(2*A))); T.setElement(5,6,(-b3/(2*A))); T.setElement(5,8,1); 
		
		return T;
	}
	
	/**
	 * Mounts the material properties matrix.
	 * @param x1 nodal coordinate of the corner 1 in direction x in cartesian system. 
	 * @param x2 nodal coordinate of the corner 2 in direction x in cartesian system. 
	 * @param x3 nodal coordinate of the corner 3 in direction x in cartesian system. 
	 * @param y1 nodal coordinate of the corner 1 in direction y in cartesian system. 
	 * @param y2 nodal coordinate of the corner 2 in direction y in cartesian system. 
	 * @param y3 nodal coordinate of the corner 3 in direction y in cartesian system.
	 * @param E young's modulus.
	 * @param v Poisson's ratio.
	 * @param t thickness of the plate. 
	 * @return B matrix.
	 */
	private IMatrix getADA(double x1, double x2, double x3, double y1, double y2, double y3, double E, double v, double t){
		IMatrix ADA = new IMatrix(7,7); 
		ADA.setZero();
		
		double A = 0.5*(x2*y3+x3*y1+x1*y2-x2*y1-x1*y3-x3*y2);
		double D0 = (E*t*t*t)/(12*(1-v*v));   double D1 = v*D0;   double Dxy = D0*(1-v)/2; 
		double i1 = A*(x1*x1 + x2*x2 + x3*x3)/12;
		double i2 = A*(y1*y1 + y2*y2 + y3*y3)/12;
		double i3 = A*(x1*y1 + x2*y2 + x3*y3)/12;
		
		ADA.setElement(0,0,4*D0*A);
		ADA.setElement(1,0,0);
		ADA.setElement(2,0,4*D1*A);
		ADA.setElement(3,0,0);
		ADA.setElement(4,0,0);
		ADA.setElement(5,0,0);
		ADA.setElement(6,0,0);
		
		ADA.setElement(0,1,0);
		ADA.setElement(1,1,4*Dxy*A);
		ADA.setElement(2,1,0);
		ADA.setElement(3,1,0);
		ADA.setElement(4,1,0);
		ADA.setElement(5,1,0);
		ADA.setElement(6,1,0);
		
		ADA.setElement(0,2,4*D1*A);
		ADA.setElement(1,2,0);
		ADA.setElement(2,2,4*D0*A);
		ADA.setElement(3,2,0);
		ADA.setElement(4,2,0);
		ADA.setElement(5,2,0);
		ADA.setElement(6,2,0);

		ADA.setElement(0,3,0);
		ADA.setElement(1,3,0);
		ADA.setElement(2,3,0);
		ADA.setElement(3,3,36*D0*i1);
		ADA.setElement(4,3,12*D0*i3);
		ADA.setElement(5,3,12*D1*i1);
		ADA.setElement(6,3,36*D1*i3);
		
		ADA.setElement(0,4,0);
		ADA.setElement(1,4,0);
		ADA.setElement(2,4,0);
		ADA.setElement(3,4,12*D0*i3);
		ADA.setElement(4,4,(4*D0*i2+16*Dxy*i1));
		ADA.setElement(5,4,(4*D1+16*Dxy)*i3);
		ADA.setElement(6,4,12*D1*i2);

		ADA.setElement(0,5,0);
		ADA.setElement(1,5,0);
		ADA.setElement(2,5,0);
		ADA.setElement(3,5,12*D1*i1);
		ADA.setElement(4,5,(4*D1+16*Dxy)*i3);
		ADA.setElement(5,5,(4*D0*i1+16*Dxy*i2));
		ADA.setElement(6,5,12*D0*i3);

		ADA.setElement(0,6,0);
		ADA.setElement(1,6,0);
		ADA.setElement(2,6,0);
		ADA.setElement(3,6,36*D1*i3);
		ADA.setElement(4,6,12*D1*i2);
		ADA.setElement(5,6,12*D0*i3);
		ADA.setElement(6,6,36*D0*i2);
		
		return ADA;
	}
	
	/**
	 * Correct the stiffness matrix for INSANE standard (dz, Rx, Ry).
	 * @param K stiffness matrix.
	 * @return k2 correct stiffness matrix.
	 */
	private IMatrix adjustMatrix(IMatrix K){
		IMatrix k1 = new IMatrix(9,9);
		k1.setZero();

		IVector c1 = new IVector(9);
		IVector c2 = new IVector(9);
		IVector c3 = new IVector(9);
		IVector c4 = new IVector(9);
		IVector c5 = new IVector(9);
		IVector c6 = new IVector(9);
		IVector c7 = new IVector(9);
		IVector c8 = new IVector(9);
		IVector c9 = new IVector(9);

		IMatrix k2 = new IMatrix(9,9);
		k2.setZero();

		IVector l1 = new IVector(9);
		IVector l2 = new IVector(9);
		IVector l3 = new IVector(9);
		IVector l4 = new IVector(9);
		IVector l5 = new IVector(9);
		IVector l6 = new IVector(9);
		IVector l7 = new IVector(9);
		IVector l8 = new IVector(9);
		IVector l9 = new IVector(9);
		
		for (int i = 0; i < 9; i++) {
			c1.setElement(i,K.getElement(i,0));
			c2.setElement(i,K.getElement(i,1));
			c3.setElement(i,K.getElement(i,2));
			c4.setElement(i,K.getElement(i,3));
			c5.setElement(i,K.getElement(i,4));
			c6.setElement(i,K.getElement(i,5));
			c7.setElement(i,K.getElement(i,6));
			c8.setElement(i,K.getElement(i,7));
			c9.setElement(i,K.getElement(i,8));
		}
		
		for (int i = 0; i < 9; i++) {
			k1.setElement(i,0,c1.getElement(i));
			k1.setElement(i,1,c3.getElement(i));
			k1.setElement(i,2,c2.getElement(i));
			k1.setElement(i,3,c4.getElement(i));
			k1.setElement(i,4,c6.getElement(i));
			k1.setElement(i,5,c5.getElement(i));
			k1.setElement(i,6,c7.getElement(i));
			k1.setElement(i,7,c9.getElement(i));
			k1.setElement(i,8,c8.getElement(i));
		}
		
		for (int i = 0; i < 9; i++) {
			l1.setElement(i,k1.getElement(0,i));
			l2.setElement(i,k1.getElement(1,i));
			l3.setElement(i,k1.getElement(2,i));
			l4.setElement(i,k1.getElement(3,i));
			l5.setElement(i,k1.getElement(4,i));
			l6.setElement(i,k1.getElement(5,i));
			l7.setElement(i,k1.getElement(6,i));
			l8.setElement(i,k1.getElement(7,i));
			l9.setElement(i,k1.getElement(8,i));
		}
		
		for (int i = 0; i < 9; i++) {
			k2.setElement(0,i,l1.getElement(i));
			k2.setElement(1,i,l3.getElement(i));
			k2.setElement(2,i,l2.getElement(i));
			k2.setElement(3,i,l4.getElement(i));
			k2.setElement(4,i,l6.getElement(i));
			k2.setElement(5,i,l5.getElement(i));
			k2.setElement(6,i,l7.getElement(i));
			k2.setElement(7,i,l9.getElement(i));
			k2.setElement(8,i,l8.getElement(i));
		}

		return k2;
	}
	
	private IMatrix signal(IMatrix K){
		IMatrix k = new IMatrix(9,9);	k.setZero();
		IMatrix l = new IMatrix(9,9);	l.setIdentity();
		l.setElement(2,2,-1);	l.setElement(5,5,-1);	l.setElement(8,8,-1);
		
		IMatrix aux = new IMatrix(9,9);	aux.setZero();
		aux.mulTransposeLeft(K,l);
		k.mulTransposeLeft(aux,l);
		
		return k;
	}

	public String getUnknownDualVariableLabel() {
		// TODO Auto-generated method stub
		return null;
	}
}
